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APPROXIMATE QUADRATURE. 



BY PROF. ASAPH HALL. 

(1). The general expression for the area contained between two ordi- 
nates, the axis of x and any curve is Jydx, the integral being taken between 
the limits x = g, and x = h; to which correspond the terminal ordinates. 
Since it is always possible to change the limits of integration to and 1 by 
putting x = g + (h — g) t, so that 

( ydx = (h — g) ) ydt, 

we need consider only the integral J\ydt. Frequently however the integra- 
tion is impossible, or the function y is unknown, and knowing by observa- 
tion or otherwise certain values of the ordinates we are obliged to resort to 
approximate quadrature. The simplest method is this : divide the interval, 
0—1, into n equal parts and join by right lines the tops of the ordinates 
drawn at the points of division. We have now n trapezoids formed by the 
n + 1 ordinates, and taking the half sum of the sides for the height of a 

trapezoid and multiplying by the common base, J = -, we have the approx- 
imate value 

ydt = J.(%y + yi + y 2 + + y n _ 1 + |y n ). 

o 
This method is obvious and must be very ancient. It is a correct principle 

of approximation and when n is infinite gives the exact value of the in- 
tegral. 

(2). The next step in the quadrature would probably be this. Instead 
of joining the tops of the ordinates by right lines assume some curve joining 
them. The general parabolic curve 

y = a -f <*i* + «2* 2 + • • • • <Vs" 3 



f 
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can be made to pass through n + 1 points, there being n -f- 1 constants at 
our disposal. Let us take a quadratic parabola, which can be made to pass 
through three points. Its axis being parallel to the ordinates,and assuming 

the axis of a; tangent to the vertex, the equation of the parabola is y = — . 

Hence the area between the axis of x, the ordinate and the curve is — : the 

. P 

arbitrary constant being zero. The area between the first and third ordi- 

nates, the axis of x and the curve is 

( x + 2Jf—x i A 6^ + 12^ + 8 ^ AVx^ (x+J) 2 (x+2Jf~\ 
Zp ' 3' p S]_p p p J' 

or Area = %J(y + 4y t + y 2 ). 

This result has been obtained by supposing the axis of x tangent to the ver- 
tex of the parabola, but a similar formula is found if we move the axis of x 
downward a distance = b. The area will then be 

2bJ +y.[y -b + 4( yi -b) + y 2 -6] = y.(y + 4y t + y 2 ), 
in which the ordinates are now measured from the original axis of x. If we 
apply this formula to the groups y 2/i^2 : yiVzVi • ViVaUs : &c -> we have 

Cydt = y.[y +y„ + 2(y 2 + y A + y & + .••) + %i + Vz +•■ •)]• 
J o 

This is known as Simpson's rule, and is in common use for computing the 
approximate values of areas; also the tonnage of vessels, the capacity of 
casks, &c. In determining the solid content y is the surface of the section 
made by a plane perpendicular to the axis of x. 

(3). The question of quadrature like that of interpolation is indeter- 
minate and admits of an indefinite number of methods. Whatever the func- 
tion y may be if it is continuous we can express the general value of an 
ordinate by means of the known theorem for the expansion of a function, 
and since the ordinates are equidistant we have by the calculus of finite 
differences 

, . . x(x - 1) A2 . x(x-l)(x-2) A , , x(x-l)(x-2)(x-3) Ai 
y=y +xJy Q + \^ J Ao + - 1.2.3 ^ + 1.2.3.4 — y ° 

+ &c. 
This is the equation of a parabolic curve like that assumed in (2), since y , 
Jy a , &c. are constants. Multiply by dx and integrate from to 1. Then 

fydx = y +l 4r —j2 A * y ° + 21 J3y ° - m Jiy ° + m A * y ° ~ &c -' 

which will be the area comprised between the ordinates y and y x . In the 
same way the area comprised between y t and y 2 will be 



/. 



ydx = y x +3^! — to Ai + 57^1 — &c, 



2 " 12 ' /x ' 24 
and similar values for the other sections. Adding these integrals and ob- 
serving that 

4 + ^i + - + 4y— 1 =2/1 —2/0+2/2 —2/1 + • • • + y»— y— 1 = y»— y » 
Ao + Ai + • • • + A-i = 4y. — 4/o + &c., 

we have for the value of the integral between the extreme ordinates, after 
some slight reductions, the formula given by Laplace for computing the 
perturbations of a comet, and which is generally given in French and Eng- 
lish books. The first term is the sum of the trapezoids which is corrected 
by the following terms depending on the differences. This formula is not 
however so convenient as the one given by Gauss, see Werke Vol. Ill, 
p.3 28. 

(4). If we multiply the preceding general value of y by dx and inte- 
grate between the limits and n, where n denotes the number of intervals, 
we have 



"^1120 16^72 8/ y ° V720 60^96 36 10/ ^° 

/ n^ n' 17n»_5n« 137n»_n»\ J6 & 

^ 15040 288 T 720 64 ^ 1080 12/ y ° T ' 



If we put n = 2, and neglect the differences above J 2 , we have, since Jy 

= yi — yo, and J2 yo = 2/2 — 2 yi + y<» 

J 2/da; = g . (y + 4^ + y 2 j, 

and the application of this formula to the whole interval will give Simp- 
son's rule. Putting n = 3, we have the formula given by Newton, 

J ydx = g. [y + 2y x + 3y 3 + y z J, 

and which Cotes calls the pulcherrima et utilissima regula. 
If n = 6, we have 

^=6y + 18^ +27^ +24^ +^A + ^o+Tfo # 2/o^. 

As J 6 y is generally a small quantity, instead of -^, put ■££$ = ^ and 
then reducing by the formula Jy = y^ — y , &c. and denoting the com- 
mon interval by 6 we have the convenient and accurate formula given by 
Mr. Weddle, 



/ 



Jo 2/<&== ioL 2/o JrVt + y * + y& + 5 ( yi + y s+ y$) + ys s 

This formula is designed for six intervals, or seven ordinates, and we have 
to take the sum of the even ordinates and the middle ordinate, add to this 
five times the sum of the odd ordinates, and multiply the whole sum by- 
three tenths of the common interval. 

(5). The following elegant method was given by Roger Cotes in his 
Harmonia Mensurarurn, published in 1722, and is an elaboration of the 
method indicated by Newton. In the expression ydt substitute for y an 
entire function of t which shall agree with all the given values of y. If 

rp _ nt(nt — 1) (nt — 2) (w. — 3) (nt — n) 

nt — m 
M m = the value of T m for nt = m, 
the function 

m— x rp 

-i . !>m • -fur 
wi—0 JJJ -m 

reduces to y, m for the value nt = m, and has n -f 1 values common with y. 
As T m and M m are functions of t only we can compute the integral 



f. 



Ml dt ~ Amy 



r. 



> 
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and then shall have 
"i 
ydt = A y + A t y t + A 2 y 2 + .... + A n y n , 
o 
the values of A , A ly _4 2 , &c, serving for all cases. The values of the 

ratio T m -t-M m are 

T o _ (nt — 1) (nt — 2 ) (nt — 3) (nt — n) 

M ~ (-1) (-2) (-3) . . . . (-n) ~ ' 

T 1 _ nt(nt — 2) (w. — 3) (nt — n) 

M, = " " 1 . (-1) (-2) (l-») 

T 2 reft^t — 1) (w. — 3) (w. — w) 

Jf, 2.1.(— 1) (2 — ») "~ "' 

A T » ~- nt ( nt — 1) (^ — 2) (nt — n + 1) 

M n ~ n.(n — 1)(» — 2).... 1 

As an example we compute the value of _4 2 for n—\, or for 5 ordinates. 
In this case we have 

T 2 = 4t(4t — 1) (it — 3) (4. — 4) = 256i 4 — 512. 3 + 304£ 2 — 48*, 
Jf 2 =2.1.(-l)(-2)= + 4. 
Integrating the value of _" 2 from to 1 we have 



AA - 256 — 51? 4- 304 — 48 - 4- 8 

**» — 5 4 + 3 2 ~ + 15' 

or .4, = + &• 

Cotes has computed these coefficients for all cases to n = 10, or for eleven 

ordinates. The following are the values to n = 6. 

For n = l, A = A 1 =\, 

For 7i = 2, A^=A % =\, A x =%, 

Form = 3, A =A S = \, ^l 1 =^l 2 =t, 

For n=^4.,A =A i = -fa, A x =A i =\%, A i =f s , 

For » = 5, A = ^ 6 = ^ftr, ^ = ^ 4 = ft ,^ 2 = 4 3 - T \\, 

For n = 6, A = A 6 = ^fo, ^ = A 6 = &, A 2 =A 4 =^, A S = T %\. 

To find the approximate value of the integral by this method we have only 

to multiply the known value of the ordinate by the corresponding value of 

A and take the sum of the products. It will be noticed that we always 

have A n _ m — A m . 

(6). The preceding formulse are in such form that it is not necessary 
to regard the differences of the ordinates, but it is best in applying such 
methods to avoid large and irregular differences, and if necessary to dimin- 
ish the interval for which the ordinates are computed. To show by exam- 
ple the approximation of the formulse let it be required to compute the area 
bounded by the two ordinates drawn midway between the center and circum- 
ference of a circle, the diameter and the curve. Taking the radius as unit, 
the interval is \, and the values of the seven equidistant ordinates are, 

y Q =y 6 z= 1^/3 = 0.8660254 

y 1= y 5 = J^/8 = 0.9428090 

y 2 =y 4 = It/35 = 0.9860133 

The exact value of this area is 

2pV{\—3?)dx = Vx\/{1 — z 2 ) + sin" 1 *!* = ^ + |=0.9566115. 

By trapezoids area = 0.9539450 error = — 0.0026665, 
" Simpson's rule " =0.9565876 " =—0.0000239, 

" Weddle's rule " =0.9566084 " =—0.0000031, 

" Cotes' method " =0.9566099 " =—0.0000016. 

Under certain restrictions the preceding formulse become exact. Thus 
Simpson's rule gives the exact integral when the bounding curve is a quad- 
ratic parabola, and the exact cubature when the solid is bounded by certain 
parabolic surfaces. There is a peculiarity in the approximation by Cotes' 
method which is worthy of notice. It is that in passing from an odd number 



of ordinates to an even number the degree of the error of the quadrature 
does not change although its amount is diminished, but when we pass from 
an even to an odd number the error is diminished two degrees. It is there- 
fore better to use an odd number of ordinates. 

(7). Hitherto we have assumed that the line of abscissas over which 
the quadrature is extended is divided into equal parts, and that y, or func- 
tion of x, is expressed by a converging series and as an entire polynomial of 
the n th degree, so thet 

y = a + a t x + a^x 2 -f- a z x* -f- . . . . + a n x n . 
If the function is an entire polynomial of the » th degree, or less, the quad- 
rature is exact, but if there are other terms, such as a n+1 x n+1 , &c.,the error 
will be produced by these terms and will be the sum of the integrals 
ja n+1 x n+1 dx, &c. In 1814 Gauss published a memoir on approximate quad- 
rature, probably having been led to consider the subject from its use in 
computing the perturbations of the minor planets. The title of this memoir 
is "Methodus nova Integralium Valores per approximationem inveniendi." 
Werke, Band 3, p. 165. Gauss lays aside the restriction of an equal division 
of the line of abscissas and shows that by a suitable division of this line the 
degree of approximation can be nearly doubled. His method consists in 
choosing the abscissas in such a way that the polynomial may be of the 2n th 
degree and the error still remains zero. 

(8). Let us take a function of x which shall produce the given values 
of y when x has certain values. The form which Lagrange employs in his 
method of interpolation gives, denoting by a, a 1} a 2 , &c. the given values 
of x, 

{x — a x ){x — a 2 )(x — a 3 ) . . ■ ■ (x — a n ) 



y = y» 



+ Vi 



(a — a t ) (a — a 2 )(a — a 3 ) 
(x — a )(x — a 2 )(x — a s ) 



+ y 2 - 



K — a) K— a 2 ) K— a 8 ) 
(x — a)(x — d\){x — a s ) . 



(a 2 — a) (a 2 — a x ) (a 2 —a z ) . 



+ y». 



(x — a)(x — a-\)(x — a 2 ) . 



(a — a„) 

(x — a„) 



. . (a; — a n ) 



. (a 2 — a n ) 



(x — «„_!) 



,(a B — «„_!)' 



*(a„— a)(a n —a 1 )(a n — a 2 ) 
This function evidently gives y = y , y = y ly &c. for x = a, x = a 1} &c. 
Let X = {x — a){x — a l ) (x — a 2 ) . . . . (x — a n ) 

-" x + Cotf" -2 + . . . c„. 



= x n+l -f ex n + Cj*" 



The numerators in the value of y are 



X 



X 



X 



x — a x — a-, x — a 



-, &c., and 



the denominators are the values of the numerators when x — a t x = a ly 
&c, and if we denote them by M , M x , M 2 , &c, then 

_ X ,y, . X , X , X 

y ~" (x— a)M + (x — a x )M \ ' Vl + {x—a 2 jl£ 2 ' V<i + * ' " *(»— a n )M n ' % ' 

Since X = 0, for x = a, we have 

a»+i _|_ ca » + Cia »-i _}. c n = 0, 

&n&X = x n+1 —a n+1 + c(a: n — a n ) + c 1 (x n ~ 1 — a n - 1 )+ c„_ t (a; — a), 

and dividing by a; — a, 

-A_ = a;" + aaf -1 + aV" 2 + aV" 3 + a" 

x — a 

+ ex"- 1 + cax n ~* + ca?x n ~ z + ca n ~ x , 

+ c^*- 2 +c x oa; n - 3 + . . . (nor- 3 , 

+ c 2 cf"- 8 + . . . c 2 a"- 3 , 

+ ««-!• 

For a; = a, this gives 

_Z_ = (n + IK + nca"' 1 + (n— ljoio"- 2 + . . . + <v_i . 
a; — a 

dX 
Hence M is the value of -=— for x = a, and likewise M lt M 2 , &c., are 

the values for x = a 1? a; = a 2 , &c. We have also 



/ 



'i^ = i + ? + JsL + + «■ 

a; — a n -)- 1 n » — 1 

+ - + -^r + + ca- 1 

n n — 1 

+ -2l 1- +c,a H ~ 2 

n — 1 



r 7Z2 • • • + c " a 



+ &c + c„_ x . 

If we write these terms in the following order 

a" + ca"- 1 + ^cr- 2 + c 2 a ,l ~ & + . . . + e„„ 1 
+ J(o-i + c a"" 2 + cja"- 3 + . . . J- c„_ 2 ) 
+ K« n ~ 2 + o a"~ 3 + 0^-* + . . . + <V_ 3 ) 
+ • 



+ 



+ 

1 



n + 1' 



it is evident that the same quantity will be produced if in the product ari- 
sing from the multiplication of the series for X by the infinite series 

x -i + x x -z + ^-3 + ^4 + . . . . = log — * 

we reject all terms containing negative powers of x and write a for x. If 
therefore we put 

/i xdx 
= X', for x = a, 
x — a 

X' 
and denote the values of the function - for x — a, x — a ly z=a if &c, 

by A 0> A t , A%, &c, we shall have for the approximate value of the integral, 

J ydx = A y + A 1 y 1 + A z y z + . . . . A n y n . 
** o 
Hence we have the means of computing the coefficients A , A±, A 2 , &c, 

when the values of a; are given equal to a, a Xi a 2 , &c. 

(9). We have now to determine the values of x so that the polynomial 
which expresses the value of y may rise to the 2n th degree and the error of 
quadrature be zero at the same time. 
We have y =/(%), and let 

f(x) = (x — a) (x — a 3 ) (x — a 2 ) .... (a; — a n ). 
As j{x) is a polynomial of the 2n th degree, or less, if we denote by Q the 
quotient of this polynomial divided by <p{x) and by B the remainder we have 

f(x) = Q.<p(x) + E. 
We now assume 12 for /(a;), since for the values a, a i , a 2 , &c, p(x) is zero, 
and the error in the quadrature will be f\ Q.p(x)dx. Since Q does not ex- 
ceed the degree n — 1 in order that the error may be zero we determine f(x) 
by the conditions, 

j <p(x)dx = 0, f x<p{x)dx = 0, . . . I af- l (p[x)dx — 0. 

By means of a known formula of reduction, or by integrating by parts, the 
integral Jx m .f(x)dx can be reduced to multiple integrals of f{x). Thus 

j <p(x)dx = i f(x)dx 
I xf(x)dx = x I f(x)dx — I f(x)dx 2 
Cx 2 cp(x)dx = x 2 J<p(x)dx — 2xjf(x)dx 2 + 2 J f(x)dx B 
fx m <p(x)dx =x v, - 1 ff(x)dx — (m— l)af' ~ 2 J f(x)dx % +(m—l)(m-—2) 
^(a;)^a; 3 + ( — l)" 1 " J (m-l)(m-2) . . .l.J <p(x)dx m . 



Hence in order tha,tfxf(x)dx > fx s <p (x)dx, &c., become zero between the limits 
and 1 the integrals ff (x)dx, / 2 f(x)dx i . . . . J m f(x)dx m must vanish for the 
same limits. We have therefore to find a function such that the function 
itself and its 1st, 2nd, 3rd, and nth differentials vanish for x=0, and x=l. 
If fix be the function sought fix must have the factors a; n+1 and (x — l) n+1 , 
and inversely every function which has the factor x* +1 .(x — 1)" +1 , and be- 
sides only constant factors, fulfils the condition. Since <p(x) is of the (n-f- l) th 
order, fix — J n+1 <p{x)dx, is of the (2ti + 2) th order. If we put fix 
= x n+1 .(x — l) n+1 .M, where if is a constant, we have 

rv ; da?* 1 

Applying the rule for the differentiation of the product of two functions of 

1 



x and multiplying by M - 



(2»i+2)(2n+l)(2n) . . . (n+2)' 



n _ . +1 (n+lf „ , (n+iy.x* „_ t 

?W- X W+^' X i "].2.(2n+2)(2n+l)- a5 

. (n+iy.n\(n- iy 

1.2.3.(2n + 2) (2n + 1) (2») * + 

4- (_i\-+i (n+l)(n){n—l)(n —2) . . . .... 1 

^ v ; '(2n+2)(2n+l)(2«,)(2n— 1 )".'"'. . n+2" 
If now we assume for n the values 0, 1, 2, 3, 4, &c, and reject terms con- 
taining negative powers of a; the roots of this equation will give the values 
of the abscissas for the Gaussian method. These roots are all real and lie 
between and 1, and when there is an odd number the middle one is ^. 
The approximate value of the integral is therefore 

j ydx = A f(a) + A^fa) 4- A 2 f(a a ) -f . . . + A a <p(a n ). 

The following are the values of the roots and coefficients for n = 0, n — 1, 

ra = 2, and n = 3. 

For n = 0. a = 0.5 A = 1. 

For n = 1. a = 0.21132487 A = \ % 

a t = 0.78867513 A t = f. 

Forn = 2. a =0.11270167 4 = &> 

o x = 0.5 ^ = |, 

a 2 = 0.88729833 A 2 = &. 

For n = 3. a == 0.06943184 4 = 0.17392742, 

a t = 0.33000948 J. x = 0.32607258, 

a 2 == 0.66999052 ^i 2 = 0.32607258, 

a 3 = 0.93056816 A 3 = 0.17392742. 
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Gauss has computed these roots and coefficients for all values of n, from 
n=0, to «.=6, to sixteen places of decimals; and a calculation of the errors 
shows that the degree of approximation as compared with the Cotesian 
method is nearly doubled. Computing the area of article (6) by the Gaus- 
sian method, and putting n = 3, or using four ordinates, I find 

log f (a) = log y(a 3 ) = 9.9554751, log ?((/,{) — log f (a 2 ) = 9.9936327, 
log^n = 9.2403681 , log A x = 9 .5133143. 

Hence we have, 

A <p{a ) = 0.1569796 
Atf (o x ) = 0.3213268 
A zf («i) = 0.3213268 
Az<p{a 3 ) = 0.1569796 

Area = 0.9566128, and error== + 0.0000013. 
The Gaussian method has not, so far as I know, been brought much into 
use, and there are reasons for this in astronomical calculations, since in com- 
puting anomalies, radii vectores, &c, the interval of time being equal, the 
method of differences gives us a valuable check on the numerical work. 
But it seems to me that this method may be advantageously applied to se- 
ries of observations. Thus if one wishes to determine the barometric or 
thermometric curve for any place, and takes the time the sun is above the 
horizon for the interval of observation, and decides to make three observa- 
tions a day, for which n = 2, he will observe his instruments about one 
ninth of the interval after sunrise and before sunset, and at noon. Hence he 
will get as good a result as from five observations equally distributed over 
the interval. In the course of a year even the saving of labor in observing 
and reducing would be very great. While the need of a judicious arrange- 
ment of the observations is more apparent in meteorology, there can be no 
doubt that this question will at length be considered in astronomical obser- 
vations; since in astronomy the objects to be observed are constantly 
increasing in number, and this increase will necessitate a more careful expen- 
diture of the labor of the observer. 



EXTRA GTION OF ROOTS. 



BY ALEXANDER EVANS, ELKTON, MARYLAND. 

Let N = the number whose root is to be extracted, r = an approximate 
root, - = the degree of the root and 22=the true root. Then, approximately, 



